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We report the results of high resolution direct numerical simulations of two-dimensional Rayleigh- 
Benard convection for Rayleigh numbers up to Ra = 10 10 in order to study the influence of tempera- 
ture boundary conditions on turbulent heat transport. Specifically, we considered the extreme cases 
of fixed heat flux (where the top and bottom boundaries are poor thermal conductors) and fixed 
temperature (perfectly conducting boundaries). Both cases display identical heat transport at high 
Rayleigh numbers fitting a power law Nu « 0.138 x Ra" 285 with a scaling exponent indistinguishable 
from 2/7 = .2857. . . above Ra = 10 7 . The overall flow dynamics for both scenarios, in particular 
the time averaged temperature profiles, are also indistinguishable at the highest Rayleigh numbers. 
The findings are compared and contrasted with results of recent three-dimensional simulations. 

PACS numbers: 47.27.te,44.25.+f,02.70.Hm 



I. INTRODUCTION 

Convection refers to phenomena where spatial inho- 
mogeneities in an advected scalar field drive a fluid flow 
which in turn transports the scalar. Convection of vari- 
ous sorts (thermal, compositional, double-diffusive) plays 
a fundamental role in a wide range of geophysical, astro- 
physical and engineering applications. Transport prop- 
erties of convective flows are of utmost interest and are 
the focus of scientific efforts worldwide. Rayleigh-Benard 
convection, the buoyancy driven flow in a fluid layer 
heated from below, is one of the fundamental paradigms 
of nonlinear physics, complex dynamics and pattern for- 
mation pQ . Despite the great deal of effort that has been 
devoted to it, however, the bulk transport properties of 
turbulent Rayleigh-Benard convection still present chal- 
lenges for theory, simulation and experiment. 

The Rayleigh number Ra is a ratio of driving due to 
buoyancy resulting from a vertical temperature gradient, 
to damping due to the fluid's viscosity and thermal diffu- 
sion. The enhancement of vertical heat transport by the 
convectively driven flow is measured by the dimensionless 
Nusselt number Nu. The Prandtl number Pr = v/k, is 
the ratio of the fluid's kinematic viscosity to its thermal 
diffusivity. The goal of many experiments, simulations, 
theories and analyses is to discern the behavior of Nu as 
a function of Ra and Pr and geometric structure (often 
aspect ratio) of the domain. When the Prandtl number 
and aspect ratio are fixed, the Rayleigh number enjoys 
the status of the control parameter. Bulk transport in 
convective turbulence is an open problem: there is still 
no universally accepted theoretical expectation for what 
the asymptotic high Rayleigh number Nu-Ra relationship 
should be [2j [3] HJ 13 El E] , and the state of affairs exper- 
imentally is unresolved and even somewhat controversial 



Analysis has played and continues to play an impor- 
tant role in this problem. Mathematically reliable lim- 
its on convective transport derived from the fundamental 
model place constraints on theories invoking uncontrolled 
approximations or incorporating additional assumptions. 
The classic rigorous result for Rayleigh-Benard convec- 
tion between fixed-temperature no-slip plates is the scal- 
ing bound Nu < cRa 1 / 2 where the prefactor c is uni- 
form in Pr for fixed temperature and no-slip velocity 
boundary conditions [I6j fT7] , The Nu^Ra 1 ^ 2 scaling 
has been proposed in several theories [2j [3j [7], albeit 
with distinct Prandtl number dependences, while ex- 
perimental results have suggested that the asymptotic 
high-Ra exponent is somewhere between 2/7 and 1/2 de- 
pending on the Prandtl number and, perhaps, other fea- 
tures. The 1/2 scaling bound remains the best known 
rigorous estimate for the arbitrary-Pr problem, although 
Nu < (Ra x InRa) 1 / 3 holds for sufficiently high Pr at 
fixed Ra |T5] and for infinite Pr uniformly in Ra [T9] . 

One proposed explanation for discrepancies among ex- 
perimental results is the effect of finite conductivity of the 
upper and lower boundaries [20j [21] . Fixed-temperature 
(Dirichlet) boundary conditions model plates of infinite 
thermal conductivity while the limit of poorly conduct- 
ing boundaries corresponds to fixed heat flux (Neumann) 
boundary conditions. Intermediate situations are mod- 
eled by interpolating (radiation) boundary conditions. It 
is well known that the thermal boundary conditions have 
a significant effect near the convective transition, decreas- 
ing the critical Rayleigh number for the onset of con- 
vection and shifting the instability to larger scales [22] . 
Indeed, the long-wave nature of the fixed-flux linear in- 
stability has been proposed as an important effect for 
pattern selection in high-Ra turbulent convection in large 
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T* = T cold (0) or &r/dz=-0(-l) equations are 



T t + (u ■ V)T = (PrR)" 1/2 AT, (1) 

d t oj + (u • V)w = (Pr/R) 1/2 Aw - dT/dx, (2) 
A^> = w, VU=o,i = 0, 0^/as| z=o ,i = 0, (3) 



U(i) 



T* = T hot (1) or 0T*/dz = -p (-1) 

FIG. 1: 2D convection cell, periodic in the horizontal direc- 
tion, indicating the different temperature boundary condi- 
tions considered: fixed temperature and fixed heat flux (non- 
dimensional values appear in parentheses). 



aspect ratio systems [23] . 

Recent direct numerical simulations in a three- 
dimensional cylindrical cell (Pr = .7) indicated that Nu is 
supressed for finite conductivity plates above Ra = 10 9 , 
and it was suggested that typical experimental conditions 
with finite conductivity plates are closer to the fixed flux, 
rather than the fixed temperature case [24]. Rigorous 
heat transport bounds for the fixed flux case were previ- 
ously derived [25], but they scale ^Ra 1 ^ 2 as in the fixed 
temperature case, well above the simulation results for 
either boundary conditions. 

In this letter we report the results of a high-resolution 
computational study of the difference between fixed flux 
and fixed temperature Rayleigh-Benard convection in 
two spatial dimensions. We restrict attention to two 
dimensions in order to guarantee full resolution of the 
boundary layers at high Rayleigh numbers [26 . Not un- 
expectedly, the critical Rayleigh number at convective 
onset and Nu immediately above onset differ, with the 
fixed flux Nusselt number exceeding the fixed temper- 
ature Nusselt number as anticipated by linear stability 
analysis [22], but the Nusselt numbers and other quali- 
tative and quantitative features of the flows are observed 
to coincide to a very high degree of accuracy at high 
Rayleigh numbers in distinction from the computational 
results reported for three dimensions. 



II. MATHEMATICAL MODEL AND 
NUMERICAL SCHEME 

The flow is modeled by the Boussinessq approximation 
for a unit density incompressible fluid in a horizontally 
periodic two-dimensional domain of height h with no- 
slip velocity boundary conditions at the top and bottom 
plates (Fig. [I]). After a change of variables depending 
on the temperature boundary conditions, the governing 



where T the temperature, u = iu + kw = id z ip — k.d x ^ is 
the velocity field, and uo = d z u — d x w is the vorticty. The 
control parameter R in the measure of the imposed ther- 
mal forcing whose definition depends on the boundary 
conditions: 

For the case of fixed temperature at each boundary 
the space, time, and temperature scales are, respectively, 
h, {h/ag5T*) 1/2 \ and 5T* = T hot - T coid , the dimen- 
sional temperature drop, where g is the acceleration due 
to gravity and a is the thermal expansion coefficient of 
the fluid. Then the thermal forcing parameter R in ([!]) 
and ([2| is the usual Rayleigh number 

R = Ra = ag5T*h 3 /vK (4) 

with temperature boundary conditions T\ z= o = 1 and 
T\ z= i =0. For the case of an imposed vertical fixed heat 
flux ~ (3, the space, time, and temperature scales are, 
respectively, h, (otgfS)~ X ^ \ and (3h. Then in place of Q, 

R = R = ag/3h*/i>K, (5) 

with dT/dz\ z=0i i = -1. 

Denote the space-time average of a function /(x, t) by 
(/). Then for the fixed- temperature boundary condi- 
tions the Nusselt number is 

Nu = l + (PrRa) 1/2 (wT), (6) 

while for the fixed-flux case it is 

Nu = ^1 - (PrR) 1/2 (^T)^ . (7) 

Note that R = RaNu [211 [25]. 

The numerical scheme used for simulating ([l])-([3| is a 
Fourier- Chebyshev spectral collocation method in space 
with classical fourth order Runge-Kutta for the time step- 
ping [27]. Computation of the momentum equation ([2| 
and the kinematic constraint ([3| are decoupled through 
use of a high order local formula for the vorticity at 
the boundary, derived from the Neumann boundary con- 
dition dip I 'dz 1 2=0,1 = for the stream function. The 
Dirichlet boundary condition i/j\ z =o i = is imposed in 
the solution of the elliptic system which is solved by 
the matrix-diagonalization procedure [28, 29 . To ensure 
benchmark quality simulations, the grid sizes were cho- 
sen with a minimum of eight grid points in the thermal 
boundary layer, defined as the distance from the bound- 
ary at which the extrapolation of the linear portion of 
the mean profile at the boundary equals the central mean 
temperature [24] [30] . Achieving this resolution over the 




FIG. 2: Temperature snapshot for fixed flux convection at 
Ra = 1.05 x 10 10 (R = 1.07 x 10 12 ): (a) full domain, (b) zoom 
image of boundary layer to illustrate resolution. 




FIG. 3: Ra-Nu data for fixed-flux simulations. The spectral 
code was used for an aspect ratio 2 cell (Spectral-2), and the 
fourth order finite difference method [32] was used for cells of 
aspect ratio 2 (FD-2), 4 (FD-4), and 8 (FD-8). 



full range of Ra, while maintaining computational effi- 
ciency and accuracy, was made possible through the use 
of the Kosloff & Tal-Ezer mapping [31] applied to the 
Chebyshev points and a high order exponential filter to 
control aliasing errors. 



III. RESULTS AND DISCUSSION 

All results reported here are for Pr = 1. Simulations 
were first performed for the fixed heat flux case in a cell 
of aspect ratio 2. Beginning at R = 600 — below onset of 
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FIG. 4: Comparison of heat transport data for fixed- flux and 
fixed-temperature simulations in a cell of aspect ratio 2. 



convection — R was increased step-by-step by a constant 
factor, allowing at the flow to settle into a steady or 
statistically stationary dynamical state before proceed- 
ing. A snapshot of the temperature field at the highest 
Rayleigh number is shown in Fig. [ij For each value of R 
the Nusselt number Nu was measured and the effective 
Rayleigh number recovered from Ra = R/Nu. 

In order to validate the fixed flux simulations before 
moving to the fixed temperature case, the numerical ex- 
periment was repeated using a fourth order finite differ- 
ence method [32 following the same protocol for cells 
of aspect ratio 2, 4, and 8, over a more restricted range 
of Ra in order to achieve the same numerical accuracy 
as in the spectral simulations. Fig. [3] shows all these 
fixed- flux Ra-Nu data sets, noting that the simulation 
times were sufficiently long to ensure that uncertainties 
in the Nu measurements are within the size of plot sym- 
bols. For fixed flux Rayleigh-Benard convection the crit- 
ical value of R at onset depends on (decreases with) the 
aspect ratio, and the wavelength at onset is set by the 
horizontal — rather than the vertical — scale of the domain 
[22] . Hence the initial pair of convection rolls is as wide 
as the aspect ratio permitted. But at each aspect ratio, 
as the Rayleigh number was increased the roles eventu- 
ally became unstable, turbulent convection sets in, and 
the flow was observed to organize itself into pairs of tur- 
bulent aspect-ratio 2 cells like those in Fig. [2] These 
observations, together with the data in Fig. [3] suggest 
that two-dimensional fixed-flux convective turbulence is 
independent of the aspect ratio at asymptotically high 
Rayleigh numbers. 

Following the same protocol, simulations with fixed 
temperatures at the boundaries were carried out for the 
same values of Ra in an aspect ratio 2 cell. Fig. [4] shows 
the Ra-Nu data for both scenarios. It is observed that 
the data are the same, within the simulations' uncertain- 
ties, at high Rayleigh numbers. A fit of the eight highest 
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Ra data points generated with the fixed flux boundary 
conditions yields Nu « 0.138 x Ra' 285 with a scaling ex- 
ponent indistinguishable from 2/7 = .2857 .... This Nu- 
Ra relationship is in remarkably close agreement with 
the latest high-precision three-dimensional simulations 
for fixed-temperature conditions in a cylinder with in- 
sulating walls [33] which produced Nu = 0.175 x Ra" 283 
for 10 7 < Ra < 10 9 at Pr = .7. 

The quantitative correspondence between fixed-flux 
and fixed-temperature (two-dimensional) turbulent con- 
vection is not limited to the bulk heat flux. The mean 
temperature profiles are observed to converge at high 
Rayleigh numbers as well. Fig. [5] shows horizontally and 
temporally averaged temperature profiles for the two sce- 
narios. At the highest Rayleigh numbers they are ef- 
fectively indistinguishable. We have not systematically 
compared other statistical quantities, but we have noted 
striking similarities between the large scale dynamics at 
fixed flux and fixed temperature [34] . 

The Nu ~ Ra 2 ^ 7 relationship was previously observed 
with fixed temperature boundaries in two-dimensional 
simulations [35], and in three-dimensional simulations 
both with rotation [36] and without [37], albeit at sig- 
nificantly lower Rayleigh numbers. Recent direct numer- 
ical simulations at much higher Rayleigh numbers (and 
Pr = .7) reported Nu ~ Ra 1//3 over nearly four decades 
up to Ra = 10 14 in a three-dimensional cylindrical cell 
of aspect ratio 1/2 [38]. It remains to be seen whether 
the heat transport in two-dimensional Rayleigh-Benard 
convection ever deviates from the Nu ~ Ra 2 ^ 7 scaling at 
higher Rayleigh numbers. 

In summary, the high-resolution simulation results re- 
ported here suggest that the plate conductivity plays no 



significant role in the Nu ~ Ra 2//? transport law at asymp- 
totically high Ra in two dimensions with periodic side 
conditions at Pr = 1. However, it remains to be deter- 
mined how various combinations of geometry, side-wall 
conditions, plate conductivity and the spatial dimension 
affect the bulk transport in high-resolution simulations 
at asymptotically high Rayleigh numbers. 
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FIG. 5: Time averaged temperature profiles (from right to 
left): Ra= 1.03 x 10 4 , 1.03 x 10 5 , 1.21 x 10 6 , 1.22 x 10 8 , and 
1.05 x 10 10 . The discrete data indicate the spatial discretiza- 
tion used to resolve the boundary layers. 
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